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Abstract. We consider vortices in the nonlocal two-dimensional Gross-Pitaevskii 
equation with the interaction potential having the Lorentz-shaped dependence on the 
relative momentum. It is shown that in the Fourier series expansion with respect to 
the polar angle the unstable modes of the axial n-fold vortex have orbital numbers I 
satisfying < \l\ < 2\n\, similar as in the local model. Numerical simulations show 
that nonlocality slightly decreases the threshold rotation frequency above which the 
nonvortex state ceases to be the global energy minimum and decreases the frequency 
of the anomalous mode of the 1-vortex. In the case of higher axial vortices, nonlocality 
leads to the instability against splitting with creation of antivortices and gives rise 
to additional anomalous modes with higher orbital numbers. Despite new instability 
channels with creation of antivortices, for a stationary solution comprised of vortices 
and antivortices there always exist another vortex solution, comprised solely of vortices, 
with the same total vorticity but with a lower energy. 
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1. Introduction 

The Gross-Pitaevskii (GP) equation for the order parameter is derived by replacing 
the interaction potential by the Fermi zero-range pseudo-potential [TJ |2] (see also the 
review Ref. [3j. For the gaseous Bose-Einstein condensates (BECs) usage of the pseudo- 
potential is justified, since the "gaseousness parameter", ga^, where g is the particle 
density and a s is the scattering length, is of order 10~ 3 . In a cold dilute gas only binary 
collisions at low energy are relevant and these collisions are characterized by a single 
parameter, the s-wave scattering length, independently of the interaction potential. 

The predictions of the mean-field theory agree with the experimental results on BEC 
in dilute gases when both the quantum fluctuations and the cloud of non-condensed 
atoms can be neglected. In connection with gaseous BECs, the predictions of the 
GP equation were first analyzed in Ref. where insightful estimates were given and 
important scaling relations were established. The mean-field theory was subsequently 
employed for account of the dynamical properties of BEC in (i) description [H] of 
the ballistic expansion of BEC jE] and (ii) quantitative theoretical account j7j of the 
condensate excitation spectra observed in Ref. jH]. In both cases the agreement with 
the experiment was within 5% without any fitting parameters. The existence of the 
macroscopic phase in BEC, which is one of the premises for introduction of the order 
parameter, was experimentally verified in the interference picture of two BECs in a 
double- well potential 0. Importantly, the quantitative results of this experiment were 
reproduced within the GP-based mean-field theory [TO] . 

The correspondence between theory and experiment was enhanced dramatically 
with the observation of quantized vortices. Vortices are experimentally created in BEC 
using various techniques of orbital momentum transfer to the condensate, such as "phase 
imprinting" in a two-component condensate [TT] . using a laser beam stirrer [T2*l I13| ITIj. 
through the decaying solitons into vortex rings fSl US]- rotation of the magnetic 
trap [T7], and rotation of the thermal cloud during the evaporating cooling process [T5] . 
Recently, coreless vortices were produced by the phase imprinting method in a spinor 
condensate [T9j . 

The mean-field theory with account for the mechanism of the angular momentum 
transfer to the condensate describes the experiments on vortices in gaseous BECs. First 
of all the vortex state must have lower energy in the rotating frame than the condensate 
without vortices. This picture resulted in the first proposed critical rotation frequency 
fl v I2II1 for observation of vortices. However, it was shown later that the nonvortex state 
remains a local minimum at the rotation frequencies higher than Q v [21] and, hence, an 
energy barrier may prevent nucleation of vortices. Also, the vortex precession around 
the center of the trap was identified with its anomalous mode in Ref. [22J. A surface 
mode version of the Landau theory was then suggested as a mechanism for vortex 
nucleation [23] resulting in a much higher critical nucleation frequency, defined as a 
minimum over all resonances with surface modes of orbital numbers I and frequencies oof. 
Q c = mm(u>i/l). The surface mode theory is generally consistent with the experiments 
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P"H ITS] II 7j . The dynamic vortex formation theory in which rotating atomic cloud 
plays an essential role was proposed in Refs. [2HI22]- The general vortex-nonvortex state 
stability diagram for finite temperatures was obtained in Ref. j2E|- The surface mode 
spectrum at finite temperatures was computed in Ref. [27] within the Popov version of 
the Hartree-Fock-Bogoliubov theory. It was found that the thermodynamic frequency 
Q v and the critical nucleation frequency Q c of the surface-mode theory increase with the 
temperature (however, the nucleation frequency due to the quadrupole mode remains 
almost unaffected). In Refs. [2H1ES] extensive numerical simulations of the GP equation 
and the Bogoliubov linear equations in spherically and cylindrically symmetric traps 
were performed, where the critical vortex nucleation frequencies were found. 

The excitation of the surface modes is not the only possible mechanism for vortex 
nucleation. Action of a localized perturbation, a small stirrer, can allow for vortex 
formation below the nucleation threshold Q c of the surface-mode theory. In the 
experiment of Ref. a nearly pure condensate was excited by a stirrer with the 
variable characteristic size. The rapid nonresonant vortex formation was observed for 
small stirrer size at the rotation frequency lower than Q c , indicating on a local mechanism 
related to the formation of vortex- antivortex pairs (see also Ref. [BO]). In Ref. [SI] a 
classical solid-body model for the angular momentum buildup was shown to describe 
the above experiment when the density of vortices is high. On the other hand, it is 
known that in effectively two-dimensional BEC the vortex-antivortex pair formation 
is the dominant mechanism of the energy transfer from a moving object to the fluid 
at velocities above the threshold value [S21 E3]- In 3D geometry the vortex rings are 
formed instead and, depending on the drag force, the object is either slowed down or 
gets trapped inside the ring core [3i] . 

In connection of applicability of the GP equation, an important mathematical result 
was also established: the Gross-Pitaevskii functional is the asymptotic limit for gaseous 
BECs of the N-particle energy functional both in three [HH] and in two 36 j dimensions 
(see also Ref. |37j). 

The concept of quantized vortices first appeared in connection with liquid helium in 
the pioneering ideas of Onsager and Feynman f3S] (see also Refs. jSH HOI HI] ) ■ Quantized 
vortices were experimentally detected in liquid helium [12]. Similar to gases, single 
vortices and vortex arrays are observed in helium II (for instance, via the ion trapping 
technique [HU EI]), this subject is reviewed in Ref. [3l?] . 

Liquid helium, in contrast to cold alkali gases, is a dense strongly correlated system 
(the parameter ga^ is of order 0.1). In helium the interatomic distance is on the order 
of the interaction range (see, for instance, table I in Ref. |3E])- The strong dissimilarity 
between the two systems manifests itself in the vortex core size: the vortex core in gases 
is much larger than the average interatomic distance, while in liquid helium it is of the 
same order. Evidently, there is no justification for neglecting the range of interaction in 
description of the macroscopic properties of liquid helium. 

Several mean field functionals for the superfluid helium are discussed in literature 
[IHllIZllIHlliniEn] (consult also the review Ref. |35]). The simplest form of the correlation 
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energy part for liquid helium is as follows jlH] 

SM = |d 3 r|^ 2 + ^ + rf(V,) 2 }, (1) 

where g is the local density, and b, c, and 7 are the phenomenological parameters. Note 
that the last term in ([!} accounts for the nonlocal interactions in liquid helium and 
is crucial for the correct description of its macroscopic properties |4T| I49j . Instead of 
the gradient nonlocality as in (0), an ad hoc potential with fitting parameters was also 
proposed jSHj- The higher-order term g 2+1 is essential for counteraction of the non- 
physical mass concentrations in a nonlocal model with attractive-repulsive interactions 
[51 . Thus, complicated density functionals are proposed to explain the properties of 
liquid helium. However, it has been known for quite some time that some aspects 
of vortex dynamics in liquid helium, such as the annihilation of vortex rings j52j . 
the nucleation of vortices jHB], the vortex line reconnection jHl], and the superfluid 
turbulence [55] are in fact captured by the local GP model. 

In this paper we study the effect of the finite-range of interaction on vortices. Our 
motivation is that, on one hand, some properties of vortices in the strongly interacting 
systems seem to be described by the local model and, on the other hand, since vortices 
appear due to balance between the dispersion and nonlinearity, modification of the 
interaction potential may lead to different stability properties of vortices. Thus our goal 
is to find out what changes in the stability properties of vortices can be attributed solely 
to nonlocality. 

Vortices in a nonlocal model were studied recently in Refs. [SU] EH], where a 
phenomenological potential was used which brought into agreement the vortex core 
size with the healing length and allowed to reproduce the Landau phonon-roton 
dispersion curve. Our nonlocal model is obtained by substitution of the Fermi zero- 
range potential by the Yukawa finite-range repulsive potential (the Macdonald function 
in two dimensions). As we study only the effect of the interaction range, we adopt 
the simplest model of a nonlocal interaction potential without the attractive part. 
Thus, we can drop the higher-order nonlinear term used to compensate collapse due to 
attraction. To facilitate the comparison with the local theory we will keep the external 
confining potential in our nonlocal GP equation (the effects solely due to nonlocality 
are independent of the external potential). For reasons of accuracy of the numerical 
analysis, we restrict ourselves to the two-dimensional model (which corresponds to oblate 
geometry of the system). 

The nonlocal GP theory was also recently used to study the effect of nonlocality of 
the attractive atomic interactions on stability of the ground state against collapse. It 
is found that nonlocality suppresses collapse jHZ] and is responsible for appearance of 
stable self-trapped configurations [SHI IBTi] . 

The paper is organized as follows. In section [21 we argue that the simplest finite- 
range interaction potential is the Yukawa potential (in 3D) and find the corresponding 
potential in 2D. Then, in sectional we show that the unstable orbital modes of the axial 
n-fold vortex have orbital numbers I satisfying < |/| < 2\n\, while the ground state 
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is spectrally stable - precisely as in the local theory In section 0] we first numerically 
study the effect of the interaction range on the critical thermodynamic frequency Q v and 
the stabilization frequency of the 1-vortex. Then we consider the effect of nonlocality 
on the stability properties of axial 2- and 3-vortices and look for the vortex-antivortex 
solutions corresponding to the n-vortex splitting instabilities caused by the nonlocality 
for various values of the vorticity n. Section contains concluding remarks. 

2. Nonlocal Gross-Pitaevskii model 

The Fermi pseudo-potential (given by the Dirac distribution) has one parameter - the 
scattering length. Thus, the simplest model which can describe nonlocal repulsive 
interactions involves introduction of an effective potential with two parameters: the 
strength of interaction and the interaction range. The symmetry property and 
assumption of the short-range interaction allow us to single out the effective potential 
with two parameters. Indeed, the contribution to the energy functional due to the 
atomic interactions reads 

Sint = §|d 3 r| d 3 r'|M/(r)| 2 ir(r',r)|M/(r')| 2 . (2) 

Here g is the strength of interactions and the kernel is normalized as J d 3 rf^(r) = 1. 
For the spherically symmetric interaction the kernel K depends only on the relative 
distance: K = K(\r' — r|). Its Fourier image is then a real function of the squared wave 
number k: K = K(h 2 ). In the zero-range limit (the Fermi pseudo-potential) K = 1. 
In the next order of the approximation we have K — 1 — e 2 k 2 , where e is the effective 
interaction range. However, in this form, the expansion has a problem when substituted 
into the Fourier integral over all k, since the second term is unbounded (the correction 
would be the dominant term). A remedy for this is to adopt an equivalent expression 
(up to the next order in e) - the Lorentzian: 

K= (l + eV)" 1 . (3) 

The Lorentzian has been used recently as an approximation of the finite-range 
potential in Refs. jSHl EH] ■ 

Though we will use the Lorentzian in our numerical simulations, the analytical 
approach presented in this paper is valid for a certain class of interaction potentials to 
be specified in the next section. The Gaussian 

K = C e exp (-e 2 k 2 ) , (4) 

which is another common choice for the phenomenological interaction potential, belongs 
to our class. In the limit of short-range interactions all the potentials of our class can 
be approximated by the Lorentzian (JHJ). 

In three dimensions the Lorentzian corresponds to the Yukawa effective potential 
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while in two dimensions the effective interaction potential is given by 

K 2D = -^K (-) , (6) 



27T0 2 

where K (z) is the Macdonald function (see, for instance, Ref. [60j). The kernels (0) 
and © reduce to the Dirac distribution in the limit of the zero-range interaction: 

UmK 3D (r) = 6 3D (r), \m\K 2D (r) = S 2D {r). (7) 

a^O a— >0 

Both kernels (0) and © have integrable singularity at r = 0. This is evident in the 
case of the Yukawa kernel, whereas for the 2D-kernel the following asymptotics holds 
K (x) ~ ln(2/x) as x — > jHOl- Hence, due to ln(x) = o(x~ a ) for all a > as x — > 0, 
the integrability property follows. 

In the following we will also use the operator form, given as A -1 , for the effective 
potentials which then become the Green functions for the positive-definite 

operator A: 

A=l-a 2 V 2 . (8) 

The corresponding nonlocal GP equation (in the reference frame rotating with the 
frequency Q) reads 

h 2 

ihd t ^ = v 2 ^ + v^f - riL z q + ^a- 1 !^! 2 . (9) 

2m 

Here L z is the angular momentum projection on the z-axis: L z = —ih(xd y — yd x ) = 
—ihde, where 9 is the polar angle in the transverse dimensions. To facilitate the 
comparison with the local model we will use the parabolic external potential confining 
in the transverse dimensions: 

mu\ 2 V Q 2 
V = — p = ^ P ' 

where p = (x,y) and p = \p\, V = h\j±/2 is the characteristic energy of the trap, 
and the characteristic length scale is a± = y/h/ (muj±). For simplicity, we consider the 
periodic boundary conditions along the z-axis, ^(z + d) = ^(z). 

In the following it will be convenient to use the dimensionless variables: 

t=-z~t, r = — , ?/' = ^= r ^, Vl= — , L z = —, g = — 10 

Here is the number of atoms in the condensate per z-period al. Note that the 
wave function ip is normalized to 1. The interaction range is also normalized by the 
characteristic length of the transverse potential a' = a/aj_. Below we will use only the 
dimensionless form of equation (jHJ), thus we drop the primes in all the variables. 

Assuming no dependence on z, the axially symmetric n-fold vortex solution reads 
-0 = e m9 ~ 1 ' flt A(p) (here A depends on n), with the amplitude A satisfying 

L A = j-V 2 + ^ + V(p) + gF - (// + fin) | A = 0, (11) 

Fo^l-a 2 ^)" 1 ^ 2 , (12) 
here V 2 = p~ l d p pd p is the radial part of the Laplacian in 2D and 2n J °° pdpA 2 = 1. 
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3. Stability of the axial vortices in nonlocal interactions 

The stability properties of the axial n-fold vortices in the local GP equation are well 
known. In BEC with repulsive interactions the unstable orbital modes ~ e lW of the 
n- vortex in the axisymmetric trap have orbital numbers / satisfying < \l\ < 2\n\ [SI] (I 
is the angular momentum of the mode). The 1- vortex has at least one anomalous mode 
(exactly one in 2D), i.e., the mode with positive norm and negative energy, qualitatively 
predicted in Ref. [S2] and found later in Ref. [SB]. The number of anomalous modes of 
the 1-vortex depends on geometry of the trap: in a prolate trap additional anomalous 
modes are possible [221 EH] , which describe bending of the vortex lines. The detailed 
phase diagram for vortex stability in 2D was numerically obtained in Ref. |2~T) 12*6] . The 
role of the excitations bearing negative energies on vortex instability was clarified in 
Ref. [SI]. Finally, it was numerically found that there are intervals of the interaction 
strength g where the axial n- vortices with n > 2 are dynamically unstable [6*5] . 

These results concern condensates with repulsive interactions. It was also shown 
that in the attractive condensate the 1-vortex solution undergoes splitting due to 
unstable quadrupole mode [66J. 

Related issues of vortices in the Ginzburg-Landau equation (which represents 
nonconservative generalization of the GP equation) were considered in Refs. [SZl EE ISH] 
and the properties of the optical vortices are reviewed in Ref. ["T"| . 

Here we show that there is a class of interaction potentials for which the unstable 
modes of the n-fold axial vortex have orbital numbers in the interval: < |/| < 2\n\, 
i.e., exactly as in the local model. 

It is convenient to use another form of the nonlocal GP equation with an additional 
dependent variable F describing the nonlinear term (below all equations are in 2D): 

(l-a 2 V 2 )F=h/f, (13) 

{-id t - v 2 + v{p) - nL z + 9 f} v = o. (14) 

Here V = V(p) is the general confining trap which allows for the axial vortex solutions 
with the amplitude decreasing to zero as p — > oo. 

The class of interaction potentials for which the main result below is valid is defined 
as follows. The general expression for the nonlinear term F, corresponding to the general 
kernel K, reads 

r d 2 p'K(p,p'M(p')\ 2 . (15) 

Vortex solutions are possible for the kernels K = K(\p' — p\). Indeed, the quantity Fo 
(a generalization of that in (|lip) is given by the following integral 

oo ( 2tt \ 

F = j p'dp' | J de'K([p 2 + p' 2 - 2pp' cos(0' - 6)]$) 1 A 2 {p'). (16) 

I J 

Obviously, F is a function of p only and this fact allows, in principle, existence of the 
axial vortex solutions. 
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We consider the class of the interaction potentials of the type K = K{\p' — p|) for 
which the coefficients of the operator A -1 in the Fourier expansion with respect to the 
polar angle, defined as 



A- 1 / = / p'dp' 



2tt 



J d(p cos(l<p)K({p 2 + p' 2 - 2ppf cos </>]*] 



f(p'), 1 = 0,1,2,..., (17) 



i.e., the radial operators with the kernels given by the expression in the square brackets 
in (JI?j) . are all positive definite. For example, the Gaussian interaction potential 



K 



1 



exp 



P 



27ra 2 "V 2a 2 

possesses this property. Indeed, the radial operator A^ 1 is given by the expression 



Ar 1 / 



i 



p'dp' exp 



P' 2 + P 2 
2a? 



PP_ 

2a? 



f(P% 



(19) 



where h{x) is the Bessel function of the second kind. There is a sufficiently large R 
such that the sign of the quadratic form of the operator A; -1 on some / = f(p) is given 
by the sign of the finite integral 

R R 



J J dpdp'xWi (f0 x(p') 

o o 



with x = exp{— p 2 /2a 2 }pf(p). Since the Taylor expansion of the modified Bessel 
function Ii(pp' /2a 2 ) is a uniformly convergent series in powers of pp' with positive 
coefficients, the positivity of the operator A; (JTHJ) follows. 

For the interaction potential ©, corresponding to the Lorentzian, the operator A; 
is given by 



A; = 1 — a 



V 2 

p 1 

p z 



(20) 



Obviously, the radial operators A;, / = 0,1,2,..., are positive definite with respect to 
the inner product defined by the integral pdp(-). 

We assume that the system (jll)) and |T6|) with a given axisymmetric external 
potential V(p) admits axial n-vortex solutions. 

Definition We will say that the interaction potential with the kernel K(\p' — p\) satisfies 
the positivity property if the associated radial operators Aj~ if77)) . / = 0, 1, 2, 3, are 
positive definite. 

Theorem The axial n-vortex solution to the system (TH)) and (fT3j) with the interaction 
potential satisfying the positivity property can be spectrally unstable only with respect 
to perturbations which have non-zero projection on the orbital basis e lW in the interval 
< |/| < 2\n\. The ground state solution is spectrally stable. 
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Proof. Let us linearize the system (J14|) - (J15)) about the axial n- vortex solution. 
Expanding the linear corrections $ and 

$ = {A(p) + e$(p, 9, t)}e md -^\ F = F (p) + eF 1 (p, 9, t), (21) 

in the Fourier series with respect to 6, 

l>0 l>0 

inserting the result into the system (|14 p -()15 p and solving for hi yields the nonlocal 
Bogoliubov equations for (ut,vi): 

id t ( Ul )=j( Ll+ A 9 ^\ A T gAK !f lA )( Ul )=JH l ( Ul ). (23) 

Here J = diag(l, —1), the operators A^ 1 are understood as acting on the products of 
Aui or Avi, and the operators L±i are defined as 

(n + l) 2 

L ±l = -Vj + V(p) + ^-J- + gF - (p + 0(n ± I)). (24) 

p z 



The matrix operator Hi introduced by (|23|) is real and symmetric and has also the 
following property 

rH l r = H_ h r = J Y (25) 

The operator if; is related to the orbital projection of the Hessian in the rotating 
reference frame. Indeed, the solution ip = A(p)e m0 ~ lflt is the stationary point of the 
Lagrange-modified energy functional in the rotating frame 



E = S-pN 



J d 2 p[\V^\ 2 + V\ip\ 2 -n^*L g ip + ^F\^\ 2 }-fiN, (26) 

(i.e., 5E/5ip* = at the solution) where N = J d 2 p\ip\ 2 is the number of atoms. 
Equation (123(1 can be obtained from the following one 



by expansion into the Fourier series with respect to 9. Using this fact and the properties 
of the operator Hi, one can write the expansion of the Lagrange-modified energy function 
(f2T?j) in powers of e as follows 

E = E + j j d 2 p(<T , $)e- mej H ess e mej ( ^ \ + 0(e 3 ) 

= E + 27re 2 J2 /pdp«X)^ ( Ul \ + G(e 3 ). (28) 

i>o i \ Vl J 
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For stability according to Lyapunov all the terms in the series in equation (|28|) must 
be positive. The n-fold axial vortex solution is spectrally unstable if the radial linear 
problems 

1 = 0,1,2,..., (29) 

allow eigenvalues a with non-zero real part (the associated solution to the linear system 
(f2*3*j) involves the exponent multiplier e at ). 

It is helpful to remind the general result on stability of the stationary points of 
Hamiltonian systems [7Tj . It is known that (i) the eigenvalues appear in quartets 
H = {a, —cr*, —a, a*}; (ii) the second-order correction to energy taken over the subspace 
corresponding to a quartet of eigenvalues with a non-zero real part has indefinite 
signature; (Hi) the stationary point may loose the spectral stability only by collision 
of the (imaginary) eigenvalues bearing the signatures of opposite sign or by collision at 
zero cr = 0. 

To a quartet of eigenvalues {a, — cr*, —a, a*}, with a being the eigenvalue from (j2£ 
the following eigenfunctions, additional to that in ()29)1 . can be associated: 

( 3 ) ■ ( % ) • JH -' 

jh - ( % ) - ( % i • < 3 °> 

This is due to the fact that Hi is symmetric real operator possessing property (f23|) . It is 
easy to establish that all the eigenvalues from a given quartet bear the same signature. 
The energy corresponding to the eigenfrequency u> (cr = — iu) equals to uj(X\J\X) where 
(X|J|X) = 2vr/ oo pdp(|X 1 | 2 -|X 2 | 2 ). 

For the following it is convenient to get rid off the diagonal terms — fll in the 
operator JHi by expanding it as follows (here / is the unit matrix) 

J Hi = J Hi - mi. (31) 

The new operator J Hi has the same eigenfunctions, which now correspond to the shifted 
eigenvalues: b = cr — iQl. The stable orbital numbers I correspond to non- negative 
eigenvalues of J Hi J Hi (i.e., —a 2 > 0). [Here we note that the instability with a 
polynomial growth in time t is not possible in the system ([23)1. since the two eigenvalues 
a and —a coinciding at zero belong to different operators, namely J Hi and JH_i, see 
equation ()3()j1 .] We will prove the inequality in the theorem by showing that the operator 
J Hi J is positive for |/| > 2\n\, I ^ and non-negative for I = with one zero mode. 
Indeed, consider the inner product of JHiJ with X = (Xi,X 2 ) T : 

oo 

(X\JHiJ\X) = 2tt J pdp^XHU + gAAl 1 A)X 1 + X*(L^ + gAAj 1 A)X 2 
o 

- g(X*gAA^AX 2 + X^AA^AX^, (32) 
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where L-y = L±i ± Ql = L + l(l±2n)p 2 . The operator L is non-negative. The 
quickest way to see this is to rewrite L as follows (due to equation (fTTjO 

Id l2 dl 

L « = -7aT/ a T p a < 33 > 

and take into account that, in the inner product defined by the integral f pdp(-), the 
r.h.s. in (|33|) can be factorized via the integration by parts. Using non-negativity of the 
operator Lq and dropping the following non- negative term on the r.h.s. of 



oc 

J pdp {X*gAA^ 1 AX 1 + X*gAK~ 1 AX 2 - g{XlgAAj 1 AX 2 + X^gAA^AX^} 
o 

oo 

= J pdp {{XI - X* 2 )gAA^ 1 A(X 1 - X 2 )} > 0, (34) 
o 

(since A^ -1 is assumed positive in the formulation of the theorem) we arrive at the 
following inequalities 

oo 

{ X\JHiJ\X)>2,Jpdp{^\X 1 \ 2 + l l^\X 2 \ 2 } 



> (I 2 - 2\nl\)2n \ ^(\X X \ 2 + |X 2 | 2 ). (35) 



o 



[The linear space in which the inner product in inequality ()35|) exists contains the 
eigenfunctions. Indeed, (Xi,X 2 ) represents the radial part of the eigenfunction 
(Xi(p)e i( ' +n)e ,X 2 (p)e i( ^ n)e ) (see equations and (|22|l). Since the latter is regular in 
the Euclidean variables (x,y), we conclude that for non-zero n and/or I, X\ = 0(p\ n+l \) 
and X 2 = 0{p\ n - l \) as p —>■ 0, while for the ground state (n = 0) 1 = 0, hence the r.h.s. 
of (|33j) is zero.] 

Let us establish when the inner product (X\JHiJ\X) is zero. Obviously, 
(X| Jifo^l^) = for X 1 = X 2 = A. For general /, the condition X 1 = X 2 = kA 
(« is a constant multiplier) is necessary for the dropped terms X±L Xi + X 2 L X 2 and 
the r.h.s. of equation (|3*4*j) to be zero, i.e., for the first row in formula (|3*3j) to be related 
by the equality sign. But then from the r.h.s. of (|35j) (the first row) we conclude that 
/ = 0, i.e., no additional zero modes are possible. Thus, for |/| > 2\n\ and I ^ the 
operator J Hi J is strictly positive and for I = it is non-negative with one zero mode 
X = (A, A) T . Such is also the operator Hi, in which case the zero mode is (A, — A) T . 
Now the eigenvalue problem (|29|l can be rewritten as follows 






where the application of the inverse of JHiJ to the vector (X 1; X 2 ) T on the r.h.s. is 
justified since for a ^ this eigenvector is orthogonal to (A, A) T by the Fredholm 
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alternative theorem. The minimal eigenvalue of J Hi J Hi is therefore defined as 



mini — a = mm ~ , (37) 

(x\(jHijyi\x) 

where the quotient is minimized in the orthogonal complement to the direction given by 
(A, A) T . Here we remind that for the spectral stability the quotient on the r.h.s. of (|37j) 
must take only non- negative values. But the latter quotient is positive for |/| > 2\n\ and 
/ 7^ 0. For I = we have one zero mode (A, — A) T . Q.E.D. 

Remark In the course of the proof we have established positivity of the operators Hi 
(the orbital projections of the Hessian for Q — 0) with \l\ > 2\n\ and I = for the n-fold 
axial vortex solution (the zero mode of the operator Hq corresponds to the invariance 
of the system (fT3 |) -([T5 |) with respect to a phase shift). 

The established positivity property, however, does not imply the thermal or 
Lyapunov stability of the ground state for non-zero Q, since the operators Hi = Hi— Ql J 
are not, in general, positive. 

It is known that there is a threshold value of the rotation frequency above which 
the ground state loses its Lyapunov stability against the surface mode excitations by 
a perturbation breaking the rotational symmetry (see Refs. [SHIEID- Indeed, this is a 
consequence of the formula (see equation (jSU) 

(X\H t \X) = (u-ni)(X\J\X), (38) 

where uo = io~ is the eigenfrequency of the oscillations with orbital number / in the 
non-rotating system, and X is the corresponding eigenfunction. For fixed / the above 
instability threshold is given by the quotient Qi = ui/l. For a general perturbation, the 
critical rotation frequency at which the nonvortex state loses its local stability is thus 
given by the formula Q c = min(u)i/l), where the minimum is determined over the all 
resonances with the surface modes |23j . 

The Hessian corresponding to the nonvortex state possesses an additional property. 
Setting n = in equation (j^ljl we get L± t = L + l 2 /p 2 =F giving 

H_ { = Hi + 2QU. (39) 

Therefore, for the nonvortex state at zero rotation (Q = 0) the partial operators JHi 
and JH_i both have complex conjugate pairs of the eigenvalues: a = ±iou. 



4. Numerical results 



We consider the two-dimensional nonlocal GP model (|13JI - (fT4")l with the interaction 
potential corresponding to the Lorentzian. In other words, we restrict ourselves to the 
straight vortex lines, which corresponds to the oblate geometry - the simplest possible 
case. Thus, investigation of the vortex bending and formation of vortex rings is beyond 
of the scope of the present paper. Though the results below are given only for the 
parabolic external potential, we have checked that addition of a perturbation to the 
potential does not affect the conclusions on the nonlocality-induced effects (see fig 2 
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below, for instance). [Similar results on the nonlocality-induced instabilities hold also 
in a completely different external potential - the tangent- shaped trap, which models the 
weak confinement. For instance, figures 1, 3, 4, and 5 below have their counterparts for 
the tangent trap. ] 

To study n-vortex solutions and their spectra within the one-dimensional systems 
(ITTjl - ffT^j) and (|2^j) in the radial variable we have used two pseudospectral methods [72] 
based on the Fourier and Chebyshev expansions and used the Fourier-based method in 
search for the vortex configurations in 2D. For the radial variable we have used of up to 
256 grid points (which is the number of Fourier modes and the degree of the Chebyshev 
polynomial), while the 2D grid was 128 x 128. The eigenvalues and eigenfunctions of 
the nonlocal linear problem (|2*3*)) were found by using the LAPACK routines. 

To find axial n-vortex solutions and two-dimensional vortex configurations we have 
minimized the generalized Rayleigh functional (the Rayleigh functional is also known 
as the Rayleigh quotient in applications to linear systems). For the GP equation we 
define the generalized Rayleigh functional as the energy functional £(ip) evaluated at 
the normalized function: 



R&)=£(jL^ = J d 2 p{\Vf\ 2 + V\f\ 2 -nrL z f + ^F(f)\f\ 2 



/=■-*- 



(40) 



2 

(recall that our interaction coefficient g is in fact u gN" with N being a fixed number 
of atoms, see equation (fTU|) ). Here the functional F((f>) is given by equation (fTKJ) and 

The most important feature of the minimization method for linear systems based 
on the Rayleigh quotient is transferred by our definition of the Rayleigh functional to 
the nonlinear models: the numerical iterations converge only to an actual solution of 
the stationary GP equation. To show this consider the variation of the functional R{4>): 



5R 



(41) 

\4>\\ 



Now, convergence of the numerical iteration <fr n to a minimizer (p^ of the Rayleigh 
functional means, evidently, that the variation on the l.h.s in equation ()41|) tends to 
zero. To see the connection to the stationary GP equation, note that as n — > oo we have 



C n = Re / dV/, 



8f* 



- Coo (42) 

/»= 



1 4"n\\ 

(existence of this limit was confirmed numerically). Setting // = and if; = (while 
the actual solution if> is normalized to 1, the variable <fi in (14 lj) is not) we obtain the 
solution to the stationary GP equation, which can be given also as flip = 5£(ip)/5if>*, 
since in the limit n — > oo equation (jjlj) dictates 

- Oo/oo = 0. (43) 

We emphasize that in our numerical scheme the chemical potential fi is computed 
alongside with computation of the stationary solution if}, namely by the iteration: 
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C n — » /i as n — > oo. That is why only actual solutions of the stationary GP equation 
can be minimizers of the Rayleigh functional. 

Moreover, it is easy to see that the Rayleigh quotient is also the Lyapunov functional 
for the local minima (for fixed number of atoms), since the second variation of R(<p) 
is given by the second variation of 8(f) with the substitution 5f = S((p/\\<p\\). [It is 
seen that the method based on the generalized Rayleigh quotient is ideally suited for 
application to computation of the stationary solutions of nonlinear systems, it will be 
given in more detail in a separate publication [73] .] 

We have searched for the energy minimizers by starting from the topologically 
different initial functions 0o hi the minimization procedure and comparing the energy of 
the resulting stationary solutions. The numerical step <ft n — > <f) n +i was performed in the 
direction opposite to the current gradient d~R/5<p n , thus ensuring minimization of the 
energy at each step. We have observed the fast convergence of the numerical scheme 
to the stationary solutions of the GP equation (we stopped the numerical iterations 
when the norm of variation 5R/5<p and the norms of the differences <p n — <p n -i an -d 
C n — C n -i were on the order of 10~ 12 ). Finally, we have checked the accuracy of the 
grid by varying the grid size and observing the changes in the solution, in all cases the 
error of the numerical solution had negligible effect on the results. 

Let us start with the nonvortex state. One of the characteristic parameters of 
the nonvortex state is the rotation frequency Q v above which it ceases to be the 
global minimum of the energy functional (in the rotating reference frame). The 
threshold rotation frequency was numerically computed by comparing the energies of 
the nonvortex state and the 1-vortex in the rotating frame. The result is presented in 
figure 1(a) (in our dimensionless system the rotation frequency is measured in halves of 
the trap frequency, see equation (fTUj) ). It is seen that even an unjustifiably large range 
of interactions only slightly affects the threshold rotation frequency Q v . 

We found that the 1-vortex solution has only imaginary eigenvalues a = —iuj for 
all a < 1, thus it is dynamically stable in the nonlocal GP model ()13|) -(|14| ) . It turns 
out that, similar as in the local model, it has one anomalous mode (for = 0), i.e., 
the eigenvalue corresponding to the orbital operator H h with I — 1, in equation (J23*j) 
and bearing a negative energy. The frequency of the anomalous mode decreases with 
increase of the interaction range a, see figure 1(b) (here and below, we chose to show the 
modes of the operators Hi with I > 0, thus our "anomalous mode" has negative norm 
and positive frequency and is related by equation (pffij) to the true anomalous mode of 
the operator H_i). It is known that the frequency of the anomalous mode gives the 
frequency of vortex precession around the axis of the trap [22j. As it was mentioned 
above, the results on the nonlocality-induced effects do not depend on the external trap. 
To illustrate this we give the analog of figure 1(b) for the perturbed parabolic trap in 
figure 2. The 1-vortex solution suffers from a slight deformation with respect to variation 
of the interaction range as is seen from figure 3. 

We now turn to the higher axial n-fold vortices. For the local GP equation it 
is known that the intervals of the interaction strength g where the axial n-vortex is 
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Figure 1. The threshold rotation frequency above which the 1-vortex has lower energy 
(in the rotating frame) than the nonvortex state vs. the strength of interaction for 
various values of the interaction range, panel (a), and the frequency of the anomalous 
mode (for £1 = 0) of the 1-vortex vs. the interaction range for various values of the 
interaction strength, panel (b). 
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Figure 2. The frequency of the anomalous mode (for fl — 0) of the 1-vortex vs. 
the interaction range for various values of the interaction strength for the perturbed 
harmonic trap (the parabolic potential perturbed by the quartic term), the dash and 
dash-dot lines. The solid line corresponds to the harmonic trap. 
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Figure 3. The shape of the 1-vortex solution for local, a — 0, and nonlocal, a — 0.5; 1, 
interactions. The interaction strength is g — 3000. 

dynamically unstable are intertwined with the intervals of dynamical stability [65J. 
According to the main result of section EJ the axial n- vortex solution may be unstable 
only with respect to the orbital perturbations ~ e lW satisfying < I < 2n (without loss 
of generality we can consider only positive n and I). The instabilities of the n- vortices 
result in vortex splitting. This fact can be seen from the Taylor expansion with respect 
to z = x + iy and z* of the perturbed vortex solution with a perturbation having nonzero 
projection on the orbital number / (see formulae ()21j) - ()22|l ) - for < / < 2n the lowest 
power in z and/or z*, which defines the vorticity at the origin, will be \n — l\, i.e., will 
come from the perturbation. 

Variation of the interaction range a for any fixed interaction strength g (in physical 
terms, the interaction strength multiplied by the number of atoms) results in appearance 
of a finite number of instability intervals, now with respect to the interaction range, due 
to resonances of the anomalous and normal modes. We have not found resonances 
with the orbital number I = 1 (which are not forbidden). The anomalous mode with 
orbital number 1 = 1, which is present for all axial n- vortices, has the closest to zero 
frequency for any range a (Q = 0). For the 2- vortex solution the instability intervals 
vs. the interaction range are illustrated in figure 4. We found that with increase of 
the interaction strength g the instability intervals move towards smaller values of the 
interaction range. 

Nonlocality is also responsible for appearance of additional anomalous modes with 
the orbital numbers higher than the vorticity (n < I < 2n). For Q = this happens for 
the axial n- vortices with n > 2, which have more than one orbital number / in the interval 
< / < 2n, where the operators Hi introduced in ()23|) may have negative eigenvalues. 
For the 2-vortex solution the additional anomalous mode has orbital number I = 3, 
see figure 5. For the 3-vortex solution, for instance, we have observed creation of such 
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Figure 4. The real, (a) and (c), and imaginary, (b) and (d), parts of the 
eigenfrequencies corresponding to the orbital linear modes with I = 2 of the axial 
2- vortex (here Q = 0). Only the anomalous mode and the resonant normal modes are 
shown (the inserts give the signature) . The 2-vortex solution loses its stability due to 
collision resonance of the anomalous mode (solid line) and the normal modes (dashed 
lines). 



anomalous mode with 1 = 5. It should be noted that for some values of the interaction 
strength there are anomalous modes with higher orbital numbers in the local model 
as well. For instance, the axial 2-vortex solution has anomalous mode with I = 3 for 
g = 2000. For the 3- vortex solution and g = 10 4 there is anomalous mode with I = 4 
(it resonantly collides with normal modes giving rise to instabilities, see below). 

The frequency of the created anomalous mode remains close to zero, therefore no 
resonances are possible with further increase of a. The threshold value of a at which the 
zero crossing occurs decreases with increase of the interaction strength, for instance, for 
g = 10 4 the threshold is a « 0.3 as compared to a ~ 0.35 for g = 5 x 10 3 in figure 5. 

Finally, we have found that nonlocality opens new splitting channels for the axial 
n- vortex in the interval n < I < 2n. In contrast, as noted in Ref. [65J, the instability 
channels of the axial n-vortices in the local GP equation satisfy I < n with the symmetric 
splitting instability, I = n, being the strongest. The situation is different in the nonlocal 
model since there is an additional parameter - the range a. For instance, it is possible 




Figure 5. Zero crossing with creation of an anomalous mode with the orbital number 
I = 3 of the axial 2-vortex. In the picture g — 5000 and Q = 0. The insert indicates 
the avoided crossing collisions of the normal modes. 




Figure 6. The imaginary part of the eigenfrequencies corresponding to the unstable 
orbital modes with the orbital numbers I = 2, 3 and I = 4 (in the insert) of the axial 
3- vortex. Here the interaction strength is g = 10 4 . 



to find such values of the interaction strength and the range that the only dynamical 
instability of the axial n-vortex has the orbital number higher than n. This happens 
already in the case of the 3-vortex, as seen from figure 6. However, the instabilities with 
the orbital numbers I > n turn out to be weak. Similar, for the 4-vortex solution, we 
have found that for a close to 0.1 the only unstable orbital mode has the orbital number 
I = 6 with Re(w) of order 1(T 3 . 
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Figure 7. The stationary vortex solutions for a = 0.052, g = 10 4 and £1 = 0.38. 
The left panel shows the 3-vortex solution and the right one - the 5-vortex solution 
comprised of 4 vortices and 1 antivortex (in the center). 

The new channels of instability with n < I < In seem to suggest existence of new 
non-axial vortex solutions, involving combinations of vortices and antivortices, which 
could minimize the energy for some rotation frequency. However, it turns out that in 
the nonlocal GP equation, similar as in the local one, the stationary vortex solutions 
involving combinations of vortices and antivortices are never the energy minimizers 
- there always exist another stationary solution with the same total vorticity, but 
comprised of vortices only, which minimizes the energy. For instance, the instability 
with I — 4 of the axial 3-vortex solution, shown in figure 6, could result in the formation 
of a non-axial 5-vortex solution with 4 vortices and 1 antivortex. Such solution was 
indeed found by the minimization starting from the seed resembling the axial 3-vortex, 
see the right panel of figure 7. For rotation frequencies Q > 0.38 this solution has 
lower energy than the ground state, however, there is at least one another solution, the 
combination of three vortices, shown in the left panel of figure 7, which has lower energy 
then the vortex-antivortex solution for all rotation frequencies. It is interesting to note 
that both these solutions do exist in the local model (a = 0). 

More complicated combinations of vortices and antivortices were found by the 
energy minimization. However, our numerical results indicate that for a stationary 
solution involving vortices and antivortices there always exist another stationary solution 
with the same total vorticity, but comprised of vortices only, which has lower energy. 
This is true for the local and nonlocal GP equations. For instance, in figure 8 we give 
such solutions with the total vorticity 8: the left panel shows the 8-vortex solution (the 
energy minimizer) and the right one - the 16-vortex solution with 12 vortices and 4 
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Figure 8. The stationary vortex solutions for a = 0.052, g = 10 4 and ft = 0.5. 
The left panel shows the 8-vortex solution and the right one - the 16-vortex solution 
comprised of 12 vortices and 4 antivortices. 

ant ivort ices. 
5. Conclusion 

In the present paper we have made an attempt to understand what changes in the 
properties of quantized vortices can be attributed to the range of interaction, when the 
latter is not negligible. The motivation was similarity of the properties of quantized 
vortices in the gaseous BECs, where the range of interactions is much smaller than the 
size of the vortex core and, hence, negligible, and in liquid helium, where the vortex 
core size is comparable to the interaction range. 

In general, the local GP equation cannot not be justified for description of the 
quantized vortices in liquid helium. It is common approach to use a nonlocal model for 
the phenomenological description of collective phenomena in such systems. However, 
there is only one local model and infinitely many nonlocal ones, if the derivation from 
the first principles is difficult or not possible to carry out. To choose between different 
nonlocal models one is left to rely on the analysis of each term in the model to understand 
the effect of it. In this paper we have analyzed solely the effect of the finite range of 
interaction on the properties of quantized vortices. The interaction range affects the 
stability properties of the n-fold vortices creating new channels for vortex splitting. 
However it has only slight effect on the threshold frequency above which the nonvortex 
state ceases to be the global energy minimum. 

Here we note the difference in the effect of nonlocality in two-dimensional and one- 
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dimensional GP equations. In 2D there is always a threshold value of the interaction 
range for appearance of the instabilities caused solely by nonlocality, while, as it is 
suggested by numerical simulations in Ref. |74j . in ID such instabilities may set in 
beyond all orders of the range parameter. 

Finally, in spite of the fact that the finite interaction range opens new splitting 
channels with creation of antivortices, for any value of the interaction range and any 
rotation frequency for the vortex-antivortex solution there is a plain vortex solution 
of the same vorticity but with a lower energy. Thus, the annihilation of a vortex- 
antivortex pair lowers the energy, which will prevent the vortex-antivortex states from 
being observed. However, a priori, this is not evident, since there is mutual interaction 
between the individual vortices which depends on the model and is not know in detail 
at a small distance. In the local GP equation, for instance, the interaction potential is 
known only at a large distance, and it has only the logarithmic falloff according to the 
Kirchoff-Onsager formula. We can explain the insensitivity of the topological structure 
of the minimizers to the interaction range if we suppose that the leading term in the 
interaction potential is due to the phase overlap in the kinetic term (which, in fact, gives 
nothing but the Kirchoff-Onsager formula). 

In connection to this, it would be interesting to find out what changes to the 
vortex interaction will bring the inclusion of a nonlinear and/or nonlocal correction 
to the kinetic energy in the density functional. Such corrections have been suggested 
in the literature on liquid helium (see, for instance, |49j). There will be significant 
changes in the structure of the energy minimizers as is suggested by a similar but 
completely integrable model which also admits vortex solutions - the Euclidean complex 
sine-Gordon equation. In the latter model the kinetic energy term is nonlinear, which 
results in the energy degeneracy with respect to the number of vortices and antivortices 
as long as the total vorticity is fixed [73]. This is a direction for future studies. 
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